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Summary 

A Navier-Stokes equations solver based on a pressure 
correction method with a pressure -staggered mesh and 
calculations of separated three-dimensional flows are 
presented. It is shown that the velocity-pressure decoupling, 
which occurs when various pressure correction algorithms are 
used for pressure-staggered meshes, is caused by the ill- 
conditioned discrete pressure correction equation. The use of 
a partial differential equation for the incremental pressure 
eliminates the velocity-pressure decoupling mechanism by 
itself and yields accurate numerical results. Example flows 
considered are a three-dimensional lid-driven cavity flow 
and a laminar flow through a 90* bend square duct. For the 
lid-driven cavity flow, the present numerical results compare 
more favorably with the measured data than those obtained 
using a formally third order accurate quadratic upwind 
interpolation scheme. For the curved duct flow, the present 
numerical method yields a grid independent solution with a 
very small number of grid points, and the calculated velocity 
profiles are in good agreement with the measured data. 

Nomenclature 

A. coefficient for incremental velocity, (i=u, v, or w) 

D h hydraulic diameter of the curved duct 

( t , m, n) index for mesh 
n. outward normal vector, n.= (n ji ,n } 

p pressure 

Re Reynolds number 

r, radius of curvature for inner (suction side) wall 

for curved duct flow 

r a radius of curvature for outer (pressure side) wall 

for curved duct flow 

u. velocity component, u.= [u,v,w] 


Xj cartesian coordinates Xj={x,y,z] 

|i molecular viscosity of fluid 

(^,T(,C) curvilinear coordinates 

Superscripts 

n iteration count 

nb neighboring grid points; (/+l/n,n), (/-l,m,n), 

(/,m-l,n), (/,m+l,n), (/,m,n-l) and (/,m,n+l) 

* current value 

‘ incremental (or corrective) value 

Subscripts 

ij index for spacial coordinates, i~{ 1,2,3} and 

j«{ 1,2,3] 

Mathematical Symbol 
X summation 

Introduction 

In many flow problems of practical importance, the boundary 
geometries are complex, and arbitrary shaped blockages may 
exist inside the flow path. However, an original numerical 
method based on the SIMPLE algorithm (refs. 1 and 2) has 
been used to solve incompressible flows whose domain can 
be discretized using orthogonal meshes. Thus, a number of 
numerical methods to extend the pressure correction methods 
for flows with arbitrary geometries have appeared (refs. 3 to 
10). Along with the efforts, various pressure correction 
algorithms (ref. 1 1) to enhance the convergence nature itself 
and to increase the convergence rate have also been proposed. 
A finite volume method to solve the Navier-Stokes equations 
is presented in this paper. The method is based on a pressure 
correction method and is valid for both incompressible and 
compressible flows including Reynolds-averaged Navier- 
Stokes equations (ref. 12). Calculation of three-dimensional 
incompressible laminar flows using the method is presented 
in this paper. Calculations of various two-dimensional flows 
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using the present numerical method can be found in refs. 12 
to 14. Example flows considered in these references include: 
a developing laminar channel flow, a developing laminar pipe 
flow, a two-dimensional laminar flow through a 90' bend 
channel, a laminar polar cavity flow and a supersonic flow 
over a compression ramp with shock wave - turbulent 
boundary layer interaction. 

A number of grid layouts have been proposed and tested to 
identify the most suitable one to solve the Navier-Stokes 
equations defined on arbitrary, complex geometries. In Shyy 
et al. (ref. 3), the standard fully staggered mesh is used to 
solve the Navier-Stokes equations defined on curved 
geometries. This grid layout can not be used to solve flows 
inside a 90' bend duct; see Malisika and Raithby (ref. 5) for 
more details. In reference 5, the standard fully staggered 
mesh is extended in such a way that the three velocity 
components are located at all grid points except at the pressure 
grid point. In this case, the number of degrees of freedom for 
velocity is tripled while the degrees of freedom for pressure 
remains the same as that of the original fully staggered mesh. 
The accuracy of a numerical method depends not only on the 
number of velocity grid points but also on the number of 
pressure grid points. Hence the accuracy can not be improved 
as much as the tripled velocity degrees of freedom might 
suggest. A collocated mesh is used in Rhie (ref. 6), Dwyer 
and Ibrani (ref. 7), Peric et al. (ref. 8) and Majumdar (ref. 9). 
In reference 6, the velocity -pressure decoupling is prevented 
by including a fourth order artificial dissipation in the 
pressure correction equation, while in reference 7, the same 
purpose is achieved by evaluating the incremental velocities 
(u',v’,w’) at mid-sides of the control volume. In Peric et al. 
(ref. 8) and Majumdar (ref. 9), a velocity-pressure coupled 
solution is obtained by interpolating the pressure gradient 
term differently from the other terms in the discrete momentum 
equation. This method is called as the “momentum 
interpolation” scheme (ref. 9). Majumdar (ref. 9) showed that 
the momentum interpolation scheme of Peri et al. (ref. 8) may 
yield a numerical result that depends on the under- 
relaxation parameter, and proposed an improved momentum 
interpolation scheme that yields a unique solution. In 
Vanka et al. (ref. 10), the velocities are located at the same 
grid points and the pressure is located at the cetroid of the cell 
formed by the four adjacent velocity grid points (for the 
two-dimensional case). The pressure -staggered mesh for the 
three-dimensional case is shown in figure 1 for clarity in the 
following discussion. The pressure-staggered mesh was 
first used in the control volume method by Vanka et al. 
(ref. 10). They mentioned that it was not easy to obtain 
convergent solutions due to the velocity-pressure decoupling. 
The mechanism that leads to the velocity -pressure 
decoupled solution is heuristically shown in reference 5. 

In the present pressure correction algorithm, a partial 
differential equation for incremental pressure is solved to 
correct the velocity and pressure field. The mechanism that 
leads to the velocity -pressure decoupled solution in the 


pressure-staggered mesh (ref. 10) is identified in the following 
section, and it is shown that the use of the partial differential 
equation for incremental pressure eliminates the velocity- 
pressure decoupling mechanism by itself. The use of the 
same incremental pressure equation in the collocated mesh 
also yields a velocity- pressure coupled solution (ref. 15). In 
the case of an orthogonal, collocated mesh, the present pressure 
correction scheme becomes very similar to the momentum 
interpolation scheme (refs. 8 and 9), and both schemes yield 
strongly diagonally dominant systems of equations for the 
incremental pressure, see reference 15 for details. 

A few pressure correction algorithms, to enhance the 
convergence nature and to accelerate the convergence rate, 
have also been proposed and are in use to solve the Navier- 
Stokes equations. Among these are the SIMPLE-C and the 
SIMPLE-R (refill). However, it is shown in Jang etal. 
(ref. 11) that the performance of each pressure correction 
algorithm (SIMPLE-C or SIMPLE-R) depends somewhat on 
the particular flow problems to be solved. With any of these 
methods, the relative error and the mass imbalance can be 
reduced only to a few orders of magnitude (e.g., KT 4 ) 
especially when a highly graded mesh is used. In the present 
method, solving the partial differential equation for the 
incremental pressure does not require the use of any 
specialized pressure correction algorithm. It is shown in the 
“numerical results” section that the present method yields 
highly converged results (i.e., the relative error and the 
mass imbalance are reduced to order of 10*) within a 
reasonable number of iterations for the curved duct flow. It 
can also be found in reference 1 3 that the relative error can 
be reduced to the same order for a transonic flow with shock 
wave - turbulent boundary layer interaction even with the 
use of a highly stretched mesh (e.g., grid aspect ratio of 
approximately 100,000). 

Numerical Method 


The incompressible laminar flow equation are given as; 
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3xj 


- = 0 


( 1 ) 
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where the subscripts i and j denote each coordinate direction, 
and the repeated indices imply summation over the index 
unless otherwise stated. 

In the pressure correction methods, the velocities and the 
pressure are decomposed as; 

Ui = u* + u[ (3) 


2 



p = p* + p' 


(4) 


where the superscript * denotes the current velocities which 
may not satisfy the conservation of mass equation yet. The 
discrete momentum equation for u.-velocity can be written as; 

A i (/,m,n)u i (/,m,n) = X{Ajuj} --|^ + S; v , 

nb OX i 

no sum on i (5) 

where A t is the coefficient of the u.-velocity at the velocity 
grid point (I*™* 11 )* S. v is the load vector originating from the 
curvilinear grid structure, the subscript nb denotes the 
neighboring grid points, and the pressure gradient is left in 
continuous form deliberately. Substituting equations (3) and 
(4) into eq. (5) yields; 

A i ( u i * +ui ) = „b{ A i(“ i ' + “')}-^ir~ +SiV ’ 

no sum on i 

(6) 

where the grid point index (/,m,n) has been deleted in 
equation 6 for convenience. The discrete u-momentum 
equation based on the current flow variables which may not 
satisfy the conservation of mass equation can be written as; 

A[U * = z(AiUi*j--^~ + Si V , no sum on i 
nb L J ox i 


(7) 


Subtracting equation (7) from equation (6) yields; 

, 1 dp' 

Uj = — , no sum on l 

Aj 3xj 


( 8 ) 


In deriving equation (8), the summation over the neighboring 
grid points in equations (6) and (7) have been disregarded. 
Substituting equation (3) into equation (1) yields; 


du] = 3u*j 
3xj 3xj 


(9) 


Substituting equation (8) into equation (9) yields the partial 
differential equation for the incremental pressure given as; 

3 [ 1 3p' 

9xj 1 A j 3xj 


.tL 

3X; 


( 10 ) 


where the last term in equation (10) represents the mass 
imbalance. Applying the standard finite volume procedure to 
equation (10) yields the discrete incremental pressure equation. 
In this case, the value of A. at the interface of the pressure 
control volumes is obtained by averaging the Aj values defined 


at the four adjacent velocity grid points (fig. 1(c)). As all the 
central-differenced finite volume equations for diffusion 
equations are strongly diagonally dominant, the present 
discrete pressure correction equation is strongly diagonally 
dominant. The momentum interpolation scheme (ref. 8) 
may yield a numerical result that depends on the under- 
relaxation parameter (ref. 9). On the other hand, the present 
pressure correction scheme does not yield a numerical result 
that depends on the under-relaxation parameter. The pressure 
correction equation, (eq. (10)), clearly states that the 
incremental pressure is driven only by the mass imbalance. 

For clarity, the mechanism that leads to the velocity-pressure 
decoupled solution, when various pressure correction 
algorithms (refs. 1, 5 and 10) are applied to the pressure- 
staggered (or collocated) mesh, is described below. In these 
pressure correction algorithms, the influencing pressure 
nodes for each velocity grid point are identified first. In tire 
case of a pressure-staggered mesh, it is argued that the 
velocity at a velocity grid point is influenced by the eight 
adjacent pressure grid points, see figure 1(b). The pressure 
gradient at the interface of the pressure control volumes is 
obtained by a veraging the pressure gradients at the four 
adjacent velocity grid points, see figure 1(c), and that of each 
velocity grid point is given as: 


dp' = a P ' dZj 

dxj 3xi 


where, for example, 



( 12 ) 

The discrete pressure correction equation is obtained by 
substituting equations (8) and (11) into equation (9) and 
integrating it over the pressure control volume. This discrete 
pressure correction equation contains 27 entries inter- 
connecting 27 adjacent pressure grid points (for the three- 
dimensional case) and is not diagonally dominant in general. 
For uniform flow with uniform mesh, the diagonal entry of 
this pressure correction equation vanishes, and thus the mass 
imbalance for a particular pressure grid point tends to 
correct only the pressure of the adjacent pressure grid points. 
The velocity-pressure decoupling occurs in such a case, and 
a convergent solution can not be obtained. On the other hand, 
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the discrete pressure correction equation, obtained by applying 
the standard finite volume procedure to equation (10), is 
strongly diagonally dominant even for highly skewed mesh. 
For uniform flows, all the off-diagonal entries of the pressure 
correction equation vanish and only 7 diagonal entries remain 
in the present pressure correction algorithm, and this discrete 
equation is strongly diagonally dominant. The rest of the 
numerical procedures are briefly described below. 

The power law upwind differencing scheme (ref. 1) is used 
for the momentum equation. The pressure gradient in the 
momentum equation is evaluated using equation (1 1) with the 
incremental pressure replaced by the pressure. The current 
flow velocities are obtained by solving equation (10) with the 
mass imbalance evaluated from the current velocities. The 
velocity corrections are obtained using equation (8) and 
incremental pressure. The velocity and pressure are corrected 
using equations (3) and (4), and the corrected velocities and 
pressure are used to obtain new current velocities. The 
systems of discrete incremental pressure equations and the 
discrete momentum equations are solved by the Tri-Diagonal- 
Matrix-Algorithm (ref. 2). The systems of equations are 
solved iteratively until the residuals become smaller than the 
prescribed convergence criteria. Each iteration consists of 3 
sweeps for the momentum equation and 7 sweeps for the 
pressure correction equation in each coordinate direction, 
respectively. The convergence parameters used are; 


Lid-Driven Cavity Flow 

The lid-driven cavity flow is schematically shown in 
figure 2. The cavity flow at Reynolds number (based on the 
transverse velocity, v, and the length of the lid) 3200 is 
considered in this paper. The flow domain is discretized by 
48x35x35 and 78x67x67 meshes in x-, y-, and z-coordinate 
directions, respectively. The convergence criterion used is 
1.0x10"*. The residuals at the time of convergence are 
Rj-S.OxlO* 3 and R^l.OxlO"*, respectively. The calculated 
transverse velocity (v) profiles on the symmetric plane at 
y =0.075 mare compared with the measured data (ref. 16) and 
the numerical results by Freitas etal.(ref. 18) in figure 3(a). 
It is shown in the figure that the present fine mesh result 
compare more favorably with the measured data than the one 
calculated by Freitas et al. (ref. 18) using a formally third 
order accurate quadratic upwind interpolation scheme. The 
calculated vertical velocity (w) profiles on the symmetric 
plane at z-0.075 m are shown in figure 3(b). Again, the 
present fine mesh result compares more favorably with the 
measured data than that by Freitas et al. (ref. 18) as shown 
in the figure. 

Laminar Flow Through a 
90’ Bend Square Duct 


R, = X 

NC 


0U ; 


3x; 


(12a) 
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< Max< 
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for t = 1,N j 



x for k = 

(u,v,w,p) • 


(12b) 


where N c is the number of pressure control volumes; the 
subscript k denotes each flow variable; the subscript t 
denotes each grid point; N denotes the number of grid points; 
and A k n+1 denotes the maximum magnitude of the k-th flow 
variable at (n+l)th iteration. The iteration is terminated if 
either equation (12(a)) or equation (12(b)) becomes smaller 
than the prescribed convergence criterion. 


Numerical Results 

The numerical method described in the previous sections is 
tested and evaluated by solving a three-dimensional lid-driven 
cavity flow (ref. 16) and a laminar flow through a 90° bend 
square duct (ref. 17). 


The laminar flow through a 90' bend square duct is 
schematically shown in figure 4. The Reynolds number 
based on the hydraulic diameter of the duct (D h =0.04 m) and 
the bulk velocity is approximately 800. The measured data 
can be found in Humphrey et al. (ref, 17). Both grid and 
flow domain independence studies are made for this flow. 
In the first case, the upstream boundary is located at 4D h 
upstream of the curved section and the exit boundary is 
located at 10D h downstream of the curved section. The 
symmetric half of the flow domain is discretized by 
68x18x33 and 110x38x61 meshes in x-, y-, and z-coordinate 
directions, respectively. In the second case, the upstream 
boundary is located at 10D b upstream of the curved section, 
and the symmetric half of the flow domain is discretized 
by a 92x18x33 grid points in x-, y-, and z-coordinate 
directions, respectively. In each case, the velocity profile of 
a fully developed square duct flow is prescribed at the inlet 
boundary and the vanishing gradient boundary condition is 
used at the exit boundary. The numerical results obtained 
using the smaller flow domain are presented below. 

The convergence history for the 110x38x61 mesh is 
shown in figure 5. It is shown in the figure that the converged 
solution is obtained in approximately 800 iterations. In the 
finite volume methods, the discrete system of equations is 
derived by integrating the governing differential equations 
over the control volume (ref. 1). For flows with arbitrary 
geometries, the number of interpolations to obtain flow 
variables at the cell boundaries for the present pressurc- 
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staggered mesh is as small as that of any other grid layout 
discussed previously. The enhanced convergence rate is partly 
attributed to the grid layout which requires fewer interpolations. 

The calculated longitudinal velocity profiles are shown in 
figure 6. The present numerical results obtained using the 
coarse mesh are almost the same as those obtained using the 
fine mesh which has approximately twice as many grid points 
in each coordinate direction. The difference between the two 
sets of present numerical results is less than a few percents 
in most of the flow domain, which shows that the present 
numerical method yields an almost grid independent 
solution with as small as 68x18x33 grid points in the x-, y- 
, and z-coordinate directions, respectively. It is also shown 
in the figure that the calculated secondary peaks of the 
longitudinal velocity near the inner wall (r=r.) are in good 
agreement with the measured data. It is shown later that such 
secondary peaks are formed by the secondary fluid motion. 

The calculated pressure distributions on the inner and outer 
walls at the symmetric plane are shown in figure 7, where s is 
the distance measured along the center line of the duct and 
s=0 at 0=0°. It can be seen in the figure that the present 
numerical results are in good agreement with those obtained 
by Humphrey et al. (ref. 17). The projected velocity vectors 
on a plane very close to the outer wall (£~0.005D h ) are shown 
in figure 8(a) and the pressure distribution on the same 
plane is shown in figure 8(b). The plane is located so close 
to the wall that the velocity component normal to the plane 
is by far smaller than the velocity component in the figure. It 
can be seen in these figures that the velocity vectors are 
aligned in the direction of decreasing pressure, which shows 
that the fluid motion in the region very close to the wall is 
governed mostly by the pressure distribution. Due to the 
strong adverse pressure gradient, a small reversed flow 
region is formed near the comer wall and the reversed flow 
region extends up the 0*37° as can be seen in figure 8(a). 
The projected velocity vectors and the pressure distribution 
on a plane very close to the side wall (r| = 0.005D h ) is shown 
in figure 9. Again, the velocity component normal to the 
plane is negligible, and thus the fluid motion in the near- 
wall region is mostly governed by the pressure distribution. 
The small reversed flow region near the comer wall, 
caused by the strong adverse pressure gradient, can also be 
seen in figure 9(a). 

The secondary fluid motions and the pressure distributions 
at four cross-sections of the curved duct are shown in 
figure 10. At the inlet of the curved section, figure 10(a), the 
fluid particle moves upward and recirculatory (or vortex) 
motion has not been formed yet. On the other hand, the 
numerical results in Humphrey etal. (ref. 17) show that the 
recirculatory motion already exists in the same region. To 
clarify any possible dependence of the numercial results on 
the location of the inlet boundary, the same flow field is 
calculated with the inlet boundary located at 10D h upstream of 
the curved section. The numerical results obtained using the 
larger flow domain are practically identical to the one 
shown in figure 10(a). It can be seen in these figures that the 


recirculatory motion becomes stronger as the fluid particle 
travels in the downstream direction and that the pressure 
decreases almost uniformly from the outer (pressure side) 
wall toward the inner (suction side) wall. These results indicate 
that the secondary fluid motion is more strongly governed by 
the strains caused by the streamline curvature than by the 
pressure distribution. It can also be found from these figures 
that the secondary peaks of the longitudinal velocity, sec 
figure 6, are caused by the mass carried toward the inner wall 
by the secondary motion of the fluid. 

Conclusions and Discussion 

A fmite volume method to solve the Navier- Stokes equations 
and calculations of a three-dimensional lid-driven cavity flow 
and a laminar flow through 90° bend square duct arc presented. 
A partial differential equation for the incremental pressure 
and the pressure-staggered mesh form the basis of the present 
numerical method. 

The use of various pressure correction algorithm (refs. 1, 
5 and 10) for pressure-staggered mesh yields ill-conditioned 
discrete pressure correction equation, and the velocity-pressure 
decoupling is caused by the ill-conditioned discrete pressure 
correction equation. For example, the diagonal entry of the 
discrete pressure correction equation obtained by applying 
various pressure correction algorithms (refs. 1, 5 and 10) 
vanishes for a uniform flow, with the flow domain discretized 
uniformly using the pressure -staggered mesh. In such case, 
the mass imbalance at a particular pressure node corrects only 
the pressure of the adjacent pressure grid points, and a 
converged solution can not be obtained. On the other hand, 
the use of a partial differential equation for incremental pressure 
yields a strongly diagonally dominant discrete pressure 
correction equation even when highly graded and skewed 
mesh is used to discretize the flow domain (ref. 13), and it 
also yields a unique solution. The present method can also 
be extended to solve compressible flows by including 
convective incremental pressure terms into the pressure 
correction equation (ref. 13). 

The use of the same incremental pressure equation in the 
collocated mesh yields a velocity- pressure coupled solution 
even though the resulting numerical method is not as 
strongly convergent as the present method adopting the 
pressure-staggered mesh. The pressure -staggered mesh is 
preferred over the collocated mesh for its strongly 
convergent nature, and is also preferred over the other grid 
layouts (ref. 5) for its compactness. As a remark, the present 
pressure correction scheme becomes very similar to the 
momentum interpolation scheme (refs. 8 and 9) in the case of 
an orthogonal, collocated mesh (ref. 15). Even in this case, 
the present pressure correction algorithm yields a unique 
solution without the use of a specialized interpolation scheme 
as the one proposed by Majumdar (ref. 9). 

In the region very close to the bottom wall of the lid-driven 
cavity, the transverse velocity (v) obtained using the 
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quadratic upwind interpolation scheme (ref. 18) compares 
more favorably with the measured data than the present 
numerical result. However, in most of the flow region, the 
present numerical results compare more favorably with the 
measured data than the other numerical results (ref. 18). A 
similar situation can be found in a lid- driven polar cavity 
flow for which the numerical results obtained using the 
present method compare more favorably with the measured 
data than those obtained using a second order differencing 
scheme (ref. 12). These observations indicate that correctly 
resolving the pressure field is as important as, or more important 
than, using a higher order differencing scheme. 

It is shown that the present numerical method yields a grid 
independent solution for the curved duct flow with a very 
small number of grid points. It is also shown that the method 
yields a rapidly convergent solution for the fine mesh case. 
The mass imbalance decreases approximately eight orders 
of magnitude within a reasonable number of iterations for 
the fine mesh case. These errors are a few orders of 
magnitude smaller than those obtainable using various 
other pressure correction algorithms (ref. 11), especially when 
a fine mesh is used to discretize the flow domain. Such an 
improved convergence nature is also attributed to the present 
pressure correction algorithm. The domain independence 
study shows that the numerical results obtained with the inlet 
boundary located at 10D h upstream of the curved section are 
practically identical to those obtained with the inlet boundary 
located at 4 D h upstream of the curved section. Both numerical 
results show that the recirculatory motion has not been formed 
at the inlet of the curved section, while the numerical results 
of Humphrey et al. (ref. 17) shows such a recirculatory 
motion at the same location. At further downstream locations, 
the present numerical method predicts a stronger secondary 
flow motion than that of Humphrey et al. (ref. 17) does. The 
present numerical results are in very good agreement with 
the measured data and indicate that the secondary peak of 
the longitudinal velocity is formed by the mass carried toward 
the inner wall by the secondary motion of the fluid. 
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(a) 2-D projection of pressure-staggered mesh. 
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Figure 1.— Pressure staggered mesh. 
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Figure 2. — Lid driven cavity flow. 
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Figure 9. — Velocity vector and pressure contour on (£, Q plane at Figure 10. — Secondary velocity vector and pressure contour. 
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Figure 10. — Concluded. 
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